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Abstract 



We obtain an exact solution of the time-dependent Schrodinger equation for 
a two-electron system confined to a plane by an isotropic parabolic potential 
whose curvature is periodically modulated in time. From this solution we 
compute the exact time-dependent exchange correlation potential v xc which 
enters the Kohn-Sham equation of time-dependent density functional theory. 
Our exact result provides a benchmark against which various approximate 
forms for v xc can be compared. Finally v xc is separated in an adiabatic and a 
pure dynamical part and it is shown that, for the particular system studied, 
the dynamical part is negligible. 
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I. INTRODUCTION 



The Time Dependent Density Functional Theory (TDFT) maps an interacting 

time-dependent N-electron system, described by a Hamiltonian of the form 

H = E $- + E Vfr - r,) + E v ext (r h t) (1) 
i zm i<j i 

with pi = —ihVi the momentum operator of the single particle, V(rj — rj) the two-particle 
interaction potential {V(ri — r^-) = e 2 /|rj — rj\ for Coulomb interaction) and v ex t{ri,t) the 
time-dependent external potential, to a non interacting time dependent N-electron system 
having the same density p(r,t). In this formalism the new Hamiltonian, also known as the 
"Kohn-Sham" (KS) Hamiltonian, can be written as: 

H K s = J2 h Ks(r u Pi,t) (2) 

i 

where 

p? 

h KS (ri,pi,t) = — + v ext (ri,t) +v H {ri,t) + v xc ([p(r, t)]; r h t) (3) 

is the effective one-particle Hamiltonian. Apart from the external (v ext (ri, £)) and the Hartree 
(vH(^i,t) = J dr'p(r',t)/\ri — r'|) part, the potential contains an "exchange-correlation" 
(xc) term (v xc ([p(r,t)];ri,t)) that is an unknown functional of the density. In the TDFT 
formalism the wave function of the effective noninteracting system is a Slater determinant 
of iV one-particle orbitals <fi(r,t) which satisfies the equation: 

d 

h KS (r, p, t)(pi(r, t) = ih—(pi(r, t). (4) 
The particle density can then be written as: 

^MHEI^M)! 2 . (5) 

i 

As in the time independent DFT, the main problem in TDFT is to find a good approxi- 
mation for v xc ([p(r, t)]; r, t). Among the most used approximations we mention the Adiabatic 
Local Density Approximation (ALDA) |3[], which is a direct extension of the static LDA to 
the time dependent problem, and the Optimized Effective Potential approximation (OEP) 
in which v xc ([p(r, £)]; r, t) is written as a functional of the single-particle orbitals and (usu- 
ally) only the exchange part is considered. Both approximations determine v xc ([p(r, t)]; r,t) 
at time t as a function of the density (or single particle orbitals) at the same time. At- 
tempts to include the "memory" of the xc potential, i.e., its dependence on the density at 
earlier times, have been hampered by the fact that such a retarded potential is a severely 
nonlocal functional of the density, i.e., it does not possess a gradient expansion in terms of 
the density 0,0). For example an early attempt by Gross and Kohn (GK) || to incorporate 
retardation within the frame of the LDA was found to be plagued by inconsistencies, such 
as the failure to satisfy the "harmonic potential theorem" |5| and other exact symmetries 
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T0|,[7[|. Only very recently [|,0,|TTJ a consistent local approximation including retardation 



has been formulated within the frame of the current-density functional theory (CDFT), in 
which the current-density, rather than the density, is used as the basic variable. 

In practice, it is not always easy to decide which of the above approximations works 
best in a concrete application. A comparative study of the performance of different approx- 
imations in a simple and well controlled situation would be very useful. As a first step in 
this direction, we present, in this paper, an exact calculation of the xc potential for what is 
probably the simplest nontrivial model of interacting electrons in a time-dependent external 
potential. This model consists of two electrons, in two dimensions, subjected to a parabolic 
potential, whose curvature (which is always positive) is periodically modulated in time. A 
concrete realization of the model could be two electrons in a quantum dot [fT2| , [T3|| with a 
time-dependent parabolic confinement potential. We shall show that (i) the time dependent 
Schrodinger equation for this system is exactly solvable by a combination of numerical and 
analytical methods and (ii) the knowledge of the exact solution can be used to compute the 
exact xc potential. Our solution for v xc ([p(r,t)];r,t) turns out to be the time- dependent 
generalization of similar calculations recently performed in the static case |14|]. The value 
of these results lies in the fact that they provide a rigorous benchmark, against which the 
merits or demerits of various approximate theories can be assessed. 

This paper is organized as follows. In Section II we discuss the model and the exact 
solution of the corresponding time-dependent Schrodinger equation. In Section III we con- 
struct the exact xc field E xc ([p(r, £)]; r, t) = — Vv xc ([p(r, t)]; r, t) both in the TDFT and in 
the time-dependent CDFT, and we discuss the difference between the two forms. We also 
compare our result with the known static limit [I3|. In Section IV, we draw a comparison 
between our exact results and the ALDA, OEP, GK approximations as well as the new ap- 
proximation presented in ref. [|TTj] (VUC). In Section V we introduce a separation between 
the adiabatic and the truly dynamic part of E xc . We conclude with discussion and summary 
in Section VI. 



II. THE MODEL 

We consider two interacting electrons of effective mass m* in a 2-D harmonic potential 
with frequency u(t) periodic in time. The background dielectric constant is e. The cor- 
responding time dependent Schrodinger equation in atomic units (h = e/y/e = m* = 1) 
is: 

h^(V? + Vl) + \oo 2 {t){rl + r\) + i-]vfr( ri) r 2 ) = i JU( ri , r 2 ) (6) 

where ri and r 2 are the electronic coordinates and r%2 is the distance between the electrons. 
Introducing the Center of Mass (CM) and Relative Motion (RM) coordinates TZ = (ri+r 2 )/2 
and r = r x — r 2 the eq. (|J) decouples in the two equations: 

(-|V^ + Lu 2 (t)n 2 )y C M(n, t) = i^ OM (n, t) (7) 

(-Vl + ]oo 2 (t)r 2 + -)* RM (r,t) = i^ RM (r,t) (8) 
4 r at 
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where 



(r ls r 2 ) = V C m(K, t)V RM (r, t) (9) 

is the orbital part of the wave function. The spin state can be either a singlet or a triplet 
(we assume 3D isotropy for the spin S). The RM wave function must be even or odd under 
inversion r — > —r, depending on whether S = or S = 1 respectively. 

For simplicity of notation, we are using "r" to indicate the RM coordinate going back to 
"r 12 " only where needed to avoid confusion. 



A. Solution in the CM channel 
(i) General analytical solution 

The problem of a quantum harmonic oscillator with a time dependent frequency has been 
studied by several authors |16|]. Equation (^) is analytically solvable for a general (periodic 
or not) u(t). The angular momentum is a constant of motion and this allows the separation 
of angular and radial coordinates. So we obtain the radial equation: 

. Id 2 lid 9 , . o lm 2 N . . d , m . , . 

'"Jaw " TnM. + w m + IW^-^ = ! S""" (7J ' 4) ' (10) 

where 

^CmCR, t) = Xn,m(K, t)e m {&) (11) 

with 

m (0) = ^=e~ %m \ (12) 

where m is a positive integer denoting the (constant) angular momentum and $ is the angular 
coordinate of the center of mass. 

The general solution of eq. (]T0|) is given by: 



X„, m (R.0 = y 2 „, (n + m) , (^) 2 exp(i(2 n + m+l)(0(O)-0(i)) 

2-TT-exp ((-^ + i^JR'J^d^) (13) 

where X(t) is a complex solution of the classical equation of motion 

X(t) = -uj 2 (t)X(t), (14) 
X(t) = \X(t)\e i ^ t) (15) 

with a phase 4>(t) satisfying the condition 



The details of the derivation of (|T3|) are given in appendix A, where it is also shown that such 
a solution can always be constructed starting from two linearly independent real solutions 
of eq. ([b§). 

We stress that eq. (0) and (p~3|) provide a complete set of solutions of eq. (^) for whatever 
u(t), provided that the condition ( [To]) is satisfied. 

In the special case of an initial value problem, i.e., if the wave function is specified at 
t = as 

*CM(K, 0) = £ Cn,mXn,m{K, O)6 m (0) (17) 

n,m 

with Xn,m(Jt, 0)G m (i?) the eigenfunctions of H(0), the subsequent time evolution is given by 

*) = £ C„, m Xn, m (^, <)6 m (0) (18) 



with the initial condition for X(t) 



X(0) = -_ (19) 
/w(0) 



X(0) = i^foj(Q) (20) 
where w(0) is the frequency of the harmonic oscillator at the initial time, 
(ii) Floquet ansatz 

If the Hamiltonian is periodic in time (in this case if u(t + T) = u(t), where T is the 



period), we can look for a basis set of solutions satisfying the Floquet ansatz [|17H l9l|: 

tf(t + T) = e- ieT V(t) } (21) 

where e (real for bound states) are called Quasi- energies (QE). The QE are defined modulo 
Q = 2ir/T. This particular basis set has properties that are similar to those of the eigenstates 
of a static Hamiltonian [|18[ . 

In our calculations, we have chosen for w 2 (t) the form: 

co 2 (t) =a^(l + Acos(fit)). (22) 

To construct Floquet type solutions of eq.( |IO|) let us first of all define the Floquet solutions 
of the classical equation of motion (|TJJ) |2~T| , |2"2"| ] as the solutions Xp(t) having the property 

X F (t + T) = e iK X F (t). (23) 

There exist two solutions of this kind pljp2| corresponding to two eigenvalues e tKl ' 2 (with 
K\ = —K-i) either complex conjugate and lying on the unit circle of the complex plane or 
real and inverse to each other. In the former case the solutions Xp(t) remain bounded in 
time; in the latter, one of them increases exponentially for t — > oo, a phenomenon known 
as parametric resonance. The actual value of K as a function of A and Q can be calculated 
from eq. (0). In this way the A, Q plane can be separated in classically stable and unstable 
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regions. The border of the regions of parametric resonance are then defined by the condition: 
e iK = ±1. 



It is evident, from the form of the general solution ([13]), that Floquet solutions with real 
QE exist within the regions of classical stability. In these regions the two classical Floquet 
solutions are complex conjugate. This means that one of the two will always satisfy the 
condition (|T6| ) (see also the equivalent condition in appendix A). If we choose this particular 
solution as the one that determines the time dependence of Xn,m(Jl,t) in eq.flT3D, then the 
{Xn,m(R., t)} form a basis of Floquet wavefunctions with QE: 

K 

e n ,m = — (2n + m+ 1). (24) 
To derive eq. (^4|) we made use of the relation : 

iK = [ T —dt = [ T i- In \X\dt + % f ^-dt = iU(T) - 0(0)). (25) 
Jo X Jo dt Jo dt 

The last equality holds because Jq din \X\/dt ■ dt = since ln|X(t)| is periodic. At the 
border of the classical regions of instability e tK = ±1, the Floquet solutions Xn,m(R,t) — ► 
and the set of QE {e„ im } degenerates to {0} if e %K = 1 or to {0, tt/T} if e %K — — 1. 

In the remaining regions of classical instability solutions of the Floquet type cannot be 
constructed. These regions are centered around the values Q = 2u>o/k, k = 0,1,2... (with 
ujq the frequency of the unperturbed harmonic oscillator) J2"l~| , |2"2"fl , so that the value Q = 



is an accumulation point for the sequence. This means that if A > 0, it is not possible 
to perform the limit Q — > without entering regions of parametric resonance and so it 
is not possible to follow the evolution of a Floquet state from a finite Q down to 0. We 
remark that the occurrence of classical parametric resonance is related to a failure of the 
conventional Floquet theorem, which ensures the existence of a complete set of Floquet 
states. The reason is that our harmonic oscillator potential, being not bounded, gives rise to 
a strictly hermitian Hamiltonian and allows only Floquet states with real QE: evidently, such 
quasi-periodic states cannot describe the motion of an electron to larger and larger distance 
from the center that the resonance process would imply. Realistic bounded potentials avoid 
this problem by allowing the possibility of complex QE in which the electron can escape to 
infinity (ionization). 



B. Solution in the RM channel 



As we did for the CM channel, we separate the angular and radial coordinates in eq. @ 
and we obtain the radial equation: 

(~§-2 ~- r §- r + + \ + £)tMM) = i^n,i(r,t). (26) 

where 

* flM (r,*) = V„ ) i(r,*)e l (7?) (27) 

with 
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0,(0) = ~^=e- m (28) 

V 27T 

and Z a positive integer, even for S = and odd for £ = 1. Here •& is the angular coordinate 
for the RM channel. Eq. (Effl) cannot be solved analytically except in the following special 

cases: 

1. Time independent case The static limit (A = 0) has been well analyzed in 3-dimension 
P | and for certain values of the frequency Uq of the unperturbed harmonic oscillator 
it is possible to have a completely analytical solution also for the RM channel (see 
[0|). Similarly it is possible to construct an analytical solution in the 2-dimensional 



case (see appendix B). 

2. Weak correlation limit We define the weak correlation limit as the regime in which the 
Coulomb interaction is negligible compared to the harmonic confinement potential. 
This means that 

- -> (29) 
a 



where I = yh/2m*ujQ is the confinement length due to the harmonic potential and 
a = h 2 e/m*e 2 is the effective Bohr radius. In our units {% = e/y/e — m* — 1) this is 
equivalent to imposing ljq —>■ oo. In this regime the coulombic term becomes negligible 
and the RM problem becomes analytically solvable (see part A of this section). 

3. Strong correlation limit in the linear response approximation with respect to A In this 
limit the Coulomb interaction dominates the harmonic confinement potential. This 
means 

I , ^ 

- -> oo (30) 

a 

(so in our units uj — ► 0) and the two electrons can be shown to perform small oscilla- 
tions about the classical equilibrium position determined by the competition between 
electrostatic repulsion and harmonic confinement. Expanding the potential energy to 
the second order in the displacement from the classical equilibrium distance ro 3> I 
and neglecting corrections of order 1/tq to the kinetic energy one obtains the effective 
harmonic Hamiltonian 

Heff(t) = |U ^ 2 (t)(r - r ) 2 - vE extA (t)(r - r ) (31) 

where pf, = —d 2 /dr 2 , \i = 1/2 is the reduced mass, u) 2 {t) = 3uj 2 + ujf(t), u)f(t) is defined 
as 

uf(t) = cjgAcos(fit), (32) 

and E extj i(t) = —oof(t)ro can be viewed as an "external force" . Apart from time depen- 
dent phase factors (see Appendix C for these factors and for details of the derivation) 
the solution, for n = I = 0, takes the form: 
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\ff(r,t) = J T (^)i e " iAtio(t)(r " ro) e^^ {r " ro+:ro{t))2 (33) 

7T 4 (it 

with the solution of 

X = -u 2 {t)X{t), (34) 

<j>(t) its phase such that d<p/dt > and x (t) the solution of 

x = -E ext;1 (t) - u) 2 (t)x {t). (35) 

If we insert in eq. ( p3|) for xq(£) its linear response approximation expression 



and for X(t) the classical Floquet solution of (|51"D with d<j>/dt > 0, we obtain a Floquet 
solution for the RM problem. 

In the general case eq. (|2q) must be solved numerically, and, to construct E xc (r,t), the 
most natural choice is to consider the dynamical equivalent of the ground state, that is, for 
the RM channel, the "lowest" Floquet state ^° RM (r, t) = V'o,o( r ; t)/y/2n. We define this as 
the state that evolves continuously from the ground state of the static Hamiltonian as the 
amplitude A of the time dependent perturbation grows from zero. From now on we only 
consider I = (and correspondingly m = for the CM channel). 

In order to calculate this Floquet state we use its property of being an eigenstate of the 
one-period time-evolution operator U(T) (\&(r,t + T) = U(T)^(r, t)) with eigenvalue e~ ieT 
(see eq. fl2lD). The idea is to calculate the matrix {U(T)ij} in a suitable basis, diagonalize 
it and find its "lowest" eigenstate - the "lowest" Floquet state. The basis we choose to 
calculate {U(T)ij} is the set of eigenstates of a two dimensional harmonic oscillator with 
angular momentum equal to zero {Ri(r, 0)}. For a general instant t, {U(t)ij} are defined by 
the equation: 

M 

R j (r,t)=J2U(t) lJ R l (r,0), (37) 
i=i 

where the sum has been truncated for practical purposes. In our calculation M = 60 - a 
value that ensures a very good convergence of the lowest QE's. Inserting for each Rj(r,t) 
the expression ( |3"7| ) into eq. (f?5|), we find for U(t)ij a system of M first order differential 
equations. Integrating this system over one period, for each Rj(r,t), we obtain {U(T)ij}. 
Since the QE Ej are defined modulo Q [19] , it is not possible to establish from the value of 
the eigenvalues of U(T) the "lowest" one. To identify it we have instead used the property 
that for A —>■ 0, Sj — ► e!-, where e® is an eigenvalue of the static (A = 0) Hamiltonian |i~8| , |20 



In practice we have followed the evolution of the ground state energy e{j for increasing A's. 
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III. CONSTRUCTION AND CALCULATION OF THE EXACT 
EXCHANGE-CORRELATION POTENTIAL 



Construction in TDFT 

If we consider two electrons in a singlet state, the KS equations reduce to a single equation 
for the doubly occupied orbital ip(r,t) 

1 d 
(__V 2 + Vext (r, t) + v H {r, t) + v xc {[p{r, t)];r, t))<p(r, t) = i^¥>(r, t). (38) 

The KS orbital can be written as 

^(v,t) = \^(v,t)\e if ^ (39) 

and its modulus is related to the density by the equation: 

p(r,t) = 2|^(r,t)| 2 (40) 

while its phase is related to the KS velocity \ KS (r,t) by 

Vf(r,t)=v KS (r,t). (41) 

If we insert expression ( |39"D in eq . (|38D and we impose that t> xc ([p(r, t)}; r, t) is real, we obtain 
two equations, one from the real part of fl3"8"|): 

^V 2 lnp(r, t) + ||V lnp(r, t)\ 2 - v ext (r, t) - v H {r, t) - v xc {[p{r, t)}; r, t) 
-^|V/(r,t)| 2 -|/(r,t) = (42) 
and the second from its imaginary part 

V ■ V/(r, t) + V/(r, t) ■ V In p(r, t) + ^ lnp(r, *) = 0. (43) 

Eq.(f42"|) can be solved for v xc (r, t) (for simplicity of notation we have dropped the dependence 
of the xc potential on the density) and we find the following explicit expression for the xc 
electric field E xc (r,t) = — S/v xc (r, t): 

E xc (r, t) = -V(^V 2 mp(r, t) + -|V lnp(r, t)\ 2 ) - E ext (r, t) - E H (r, t) 

-V(-i|V/(r,t)| 2 -|/(r,t)), (44) 



where E H = — Vi># and E ext = —Vv ext . The last two terms of eq.(|44T) are peculiar of the 
time dependent problem while the first three correspond to the static expression |14| for 
E 

E £«c( r ) = _ V (Iv 2 lnp(r) + ^|Vlnp(r)| 2 ) - E ext (v) - E*(r). (45) 
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Eq.(f43|) is a first order partial differential equation for V/ and is equivalent to the continuity 
equation for the non interacting KS system. It shows that v^s(r, t) = V/(r, t) is in general 
a non trivial functional of the density. We stress that v^s(r,t) is in general not the same 
as the exact velocity field; only the longitudinal part of the KS current must coincide with 
the longitudinal part of the physical current due to the continuity equation. 
Construction in time dependent CDFT 

The time dependent CDFT differs from the TDFT in that not only the density but also 
the current density calculated from the KS single particle orbitals is exact. 

In order to accomplish this, one introduces an xc vector potential, A xc in the Kohn-Sham 
equation 0,0]. The KS Hamiltonian J7| is now: 

H KS {t) = J2(liPi + A *c(ri, t)f + v H {r h t) + v ext {r h t)) (46) 
i 1 

which yields both the correct density and current. In the case of two electrons in a singlet 
state, we get, for the occupied orbital <p(r,t), the time dependent Schrodinger equation: 

1 d 
((-(p + A xc (r, t)f + v H (r, t) + v ext {r, t)))<p(r, t) = i^(r, *) (47) 

and we can now follow the same procedure used in TDFT to find an explicit expression for 
E xc (r, t). We obtain: 

E«.(r,t) = -A xc (r,t) 

= -V(^V 2 lnp(r,t) + i|Vlnp(r,t)| 2 ) - E« rt (r,t) - E H {r,t) 

+V(^ 2 )+v, (48) 

where v(r, t) is the exact velocity of the interacting system, v(r, t) = V/(r,t) + A(r,t). 
From the imaginary part of eq. fl47l) (or equivalently from the continuity equation) we get 
the first order partial differential equation 

d 

V ■ v(r, t) + v(r, t) • Vlnp(r,t) + — lnp(r,i) = 0. (49) 

The advantage of this formulation is that it expresses E xc (r,t) as a function of the 
physical quantities v(r, t) and p(r,t). 
Circularly symmetric states 

If the time dependent state is circularly symmetric, as in the case we are studying, then 
the current is purely radial, therefore purely longitudinal, and the two expressions (44|) and 
(f48|) coincide (that is v(r,t) = VA's(r,t)). Thus, in this case, there is no difference between 
the time dependent DFT and CDFT. Eq. (f43|) can be easily integrated yielding 

9 N 1 f r dpir'.t) . , . 

^^-^mXt^' (50) 

Linear response 
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In the limit of small external time-dependent perturbation (A — > in eq. (p2|)), we 
expand all the quantities to first order in A, i.e. E exf = E ext fi + E ea;t>1 , E H = E# )0 + E H>1 , 
E xc = Ea-co + E xc l and hip = lnp + Pi/ Po where the subscripts "0" and "1" indicate 
respectively zero and first order with respect to A. Then from eq. ([44"D we obtain 

E xc , Q (r) = -V(ivinp + i|Vlnp | 2 ) - E extfi (r) - E H>0 (r) (51) 

and 

E xc>1 (r,t) = V-/(r,t) - -V(V 2 ^ + Vlnp • V^) - E ext ^r,t) - E H>1 (r,t), (52) 

at 4 po Po 

where we have neglected also the terms of order v 2 (r,t) = |V/(r, t)\ 2 . 
Calculation of the exact E xc (r, t) 

We have considered in detail two sets of the parameters ujq, Q and A that appear in 
eq . (|22"D , corresponding to high and low correlation. The two sets are given in Table |. The 
values of lvq have been chosen such that it is possible to construct analytically the solution 
of the corresponding static Schrodinger equation (see appendix B). The values of A and fl 
have been chosen so that the system is in the linear response regime but well above the 
regions of parametric resonance for the CM channel and above the first excitation energy of 
the system. 

The exact time dependent densities are plotted in Fig.l and Fig. 2. For the weak correla- 
tion parameter ujq = 1 the density is centered at the origin as we can expect from the exact 
solution in the weak correlation limit (a Gaussian centered at the origin, see section II). In 
the case of high correlation (ujq m 0.02), on the other hand, the maximum of the density is at 
finite distance from r = in agreement with the form that the RM wave function assumes in 
the strong correlation limit eq.(|33|) (an annulus of average radius r /2, with r the classical 
equilibrium distance of the two electrons): the increased strength of the Coulomb repulsion 
in respect to the harmonic confinement pushes on average the two electrons far from each 
other. In these plots the solid line represents the static limit while each of the broken lines 
corresponds to the time dependent p(r, t) at different times. 

In the insets of Fig.l and Fig.2 we plot the time dependent velocity v(r, t) = V/(r, t), 
that is the other necessary ingredient to calculate E xc (r, t). As the plots show, the motion 
is approximately a "breathing" motion: the velocity is zero at the origin while for r ^ it 
increases almost linearly. The asymptotic behavior is linear in r with a correction in 1/r 2 . 
In Fig. 3 (ujq = 1) and Fig.4 (u ~ 0.02) we finally plot the results for the field E xc (r, t). The 
solid line represents the static limit while each of the broken lines corresponds to E xc (r,t) 
at different times. We choose to plot this quantity instead of the more traditional v xc {r,t) 
since this is the meaningful physical quantity whose asymptotic behavior does not depend on 
arbitrarily fixed time dependent constants. As can be seen from the plots, as the correlation 
in the system increases, the positive peak of E xc (r, t) for small r increases too. This is related 
to the enhancement of the strength of the Coulomb repulsion that, as seen before, pushes 
the maximum of the density away from the origin. E xc (r, t) can be viewed as a force in the 
KS system which, where positive, contributes to drag the particles away from the origin. 

Starting from the asymptotic form of the RM wave function for r±2 — ► oo, that is, 
^rm oc rf 2 ■ (1 + (6sr + ibcx)/r) ■ exp(-r 2 2 0/4), with a real, 6 Q = — (6^ + 6 R 0/(20))/0 and 



11 



defined by = —u 2 (t)b^ + <ft, we get, with a little algebra, the asymptotic form of 
the density (p(r,t) oc r 2a (l + ■ exp(— r 2 <p)) and of the velocity (V/(r, i) ~ — r/2 ■ 

<i(ln0)/<it — b-^/r 2 )). From these behaviors, using eq. ([44]), we can derive the asymptotic 
behavior for the xc field that results to be the same as in the static case: E 2C m — 1/r 2 for 
r — > oo. 

In Fig. 3 and Fig.4 we also plot, for comparison, the Local Density Approximation (LDA) 
of the static E xc (r). As can be seen from the figures, LDA behaves, in general, reasonably 
well except for small r (for which in the weak correlation case has even the wrong sign) and 
for large r (for which decreases exponentially). 



IV. COMPARISON WITH APPROXIMATE THEORIES 

In this section we discuss the comparison between our exact result for ~E xc (r,t) and the 
results obtained from the most used approximations, namely the Adiabatic Local Density 
Approximation (ALDA), the Optimized Effective Potential (OEP), the Gross and Kohn 
(GK) approximation, and the hydrodynamic approximation recently introduced by Vignale, 
Ullrich, and Conti (VUC). The expressions for the xc electric field in these approximations 
are: 

• ALDA [§: 

Vi LDA (r,t) = -V(^(pMK. c (p))) (53) 

where e xc (p) is the xc energy per particle of the homogeneous electron gas. In 2D it is 
given by: 

£xc = —z + "7T ~ w ■ -T2— , (54) 

<yKT » A 1 + ai^T s + a 2 r s + a 3 rg n e 

' (55) 



'Tip 

a = -0.3568, a x = 1.1300, a 2 = 0.9052, a 3 = 0.4165 [g|,[24]] 

OEP approximation (which in this simple case is equivalent to the Hartree Fock ap- 
proximation) HQ 

E° c EP (r,t) = -V(- 1 -v H (r,t)) (56) 

GK approximation || (valid in the linear response regime): 

Eg(r,w) = -V( Pl (r,u;)f xc (p ,u;)) (57) 

where E XC) i(r, uS) is a Fourier component of E XCjl (r, t) (defined, with pi(r,t), in the 
"Linear response" section) and f xc ,L(Po,u) is the longitudinal part of the frequency 
dependent xc kernel of the homogeneous electron gas. 
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VUC approximation [TT[] for a circularly symmetric potential in 2D (valid in the linear 
response regime): 

B"F(r, U ) = E£?»(r, U ) + ^(pKUM - %)(^ + ^>)) 
2 

--^(Pofxc,T(po,u))vi(r,uj)) (58) 

where t>i(r, t) is the velocity field and f xc ,L(po,u), f X c,T(po,w) are the longitudinal and 
the transverse part of the frequency dependent xc kernel of the homogeneous electron 
gas. 

In both GK and VUC we have used for f xc (po, u) the expressions recently obtained by 



Nifosi et al. [25 



The comparison between the exact E xc (r, t) and its approximations is made plotting its first 
Fourier component E xc (r, Q). Since our calculations were done in the linear response regime, 
the difference between E xc (r, Q) and ~E XCj i(r, Q) is negligible. 

In the case of weak correlation (Fig. 5, uo = 1) the ALDA, reproduces, except for 
small r, the general trend, though underestimating the peak of the potential. In the strong 
correlation case (Fig. 6, u ~ 0.02) it underestimates the potential for very small r, but, 
for intermediate values, it gets closer to E xc (r, Q). For weak correlation and for small values 
of r the OEP approximation does not reproduce the exact behavior, while in the region in 
which E xc (r, fl) is significantly non zero it gets closer to the exact result. In the limit of 
zero correlation the OEP, which is equivalent to Hartree-Fock theory for this system, would 
give the exact result. Its behavior gets worse when the correlation increases (see Fig. 6): it 
is the only approximation that does not even reproduces the first peak of E xc (r, Q). On the 
other hand, this is the only approximation that has the correct asymptotic behavior — 1/r 2 
for r — > oo. In the weak correlation case (Fig. 5) the GK approximation has a behavior 
similar to the OEP (except for the asymptotic behavior that is not reproduced correctly), 
while, for strong correlation (Fig. 6) it reproduces the correct trend but underestimates 
E xc (r, Q) for small values of r and overestimates it for intermediate values. In the case of 
weak correlation (Fig. 5) the VUC approximation does not reproduce the exact trend for 
small r, while, for intermediate values it get closer to E xc (r, fl) though underestimating its 
peak. For strong correlation its behavior is almost indistinguishable from the ALDA. 



V. "ADIABATIC" AND "DYNAMIC" EXCHANGE CORRELATION 

POTENTIALS 

There has been considerable effort, in recent years, aimed at the construction of a fully 
dynamic xc potential, which, unlike the ALDA potential, should depend on the density 
at all previous times, i.e., have a memory. In order to assess the importance of these 
"memory effects", we shall now separate E xc (r,t) in an "adiabatic" part E^(r, t) contain- 
ing the adiabatic evolution of the static exact E xc (r) and a "dynamical 1 one E^(r, t), 
peculiar of the time dependent problem. Comparing equations (0) and (f45[) , it is easy 
to identify — V(V 2 In p(r, t)/A + |Vlnp(r, t)| 2 /8) — E H (r,t)) as adiabatic terms, while 
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— V(— |V/(r, t)| 2 /2 — <9/(r, t)/dt) are peculiar to the time dependent case. However E ext (r, t) 
is not an explicit functional of the time dependent density and must be treated more care- 
fully. The adiabatic part of E exi (r, t) is defined as the electric field E^ t (r, t) which, when 
applied to the physical interacting system, would yield the exact density p(r, t) while the 
system remains in the instantaneous ground state. 

It is then possible to define the " dynamicaT part of the external field as what remains 
after subtracting the adiabatic part: 

E*(r,t)=E ext (r,t)-E^(r;t). (59) 

Now we can separate, in the case of the two electron problem, E xc (r, t) in an adiabatic 
(Ef c (r,t)) and dynamical (E^(r,t)) part, E xc (r,t) = Ef c (r,t) +E^(r,t), where 

E*(r,t) = -V(-5£M _ i(V/M)) 2 ) - E* (r,f) (60) 

E2(r,t) = -V(^Vln(p(r,t))) 2 + iv 2 (ln(p(r, t)))]) 

-E H (r,t)-E^(r,t). (61) 
In the linear regime, using the linearized expression (|52| ) we obtain: 

KtAr,t) = -V(V^ + ^Vlnp V^) - E^r.t) - E Hfl (r,f) (62) 
4 p 4 p 

EjjCr.^sE^r.^-^r.t) (63) 

= v--E ext7l (r,t) + -E? xttl (r,t) (64) 

where we have used the fact that V/(r, i) = v(r, t) and neglected terms of order f 2 (r, t). 

The difficulty in the calculation of E^ f (r, t) is that in general its form is unknown and 
leads to a non separable, two-electron Schrodinger equation. In our case it is possible to 
calculate analytically E^ t (r, t) and its counterpart E* f (r, t) in the limit of extremely weak 
and extremely strong correlation, but for a general set of parameters it will be necessary to 
find an approximation for E"i(r, t). 

We will now show that ~E xcl {r,t), in this system, vanishes exactly in both the weak and 
strong correlation limits, and it is likely to be very small in the intermediate cases. 

Calculation of Ej^r,*) 

"Weak" correlation limit 

In this regime the response of our system to the external potential 

v ext (r,t) = -(u 2 +ul(t))r 2 (65) 

with ui\{t) defined by eq. fl32j), is a "breathing motion", i.e. it can be described as a periodic 
transformation with a length scale 

i(t)oc(^(t))-i<x|X(0|. (66) 

Then we can calculate explicitly all the quantities appearing in eq. (|64"D . The velocity 
field is given by 
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*<'.«>=j§'=<ll*MI>» < 67 > 

v( r ,t) = (^ln(|X(t)|))f (68) 
the external field is given by 

E extA {r,t) = -ul(t)r, (70) 
and the "adiabatic" external field is given by 

E^, 1 (r,t) = -(|0(t)) 2 r. (71) 

This is also the exact E^ tl (r, t) corresponding to the noninteracting harmonic oscillator 
problem. 

We can now prove that ~E x y cl (r, t) = in this approximation. Using in eq.(69) |X(t)| = 
X(t) exp(— i<f)(t)), X = —oj(t) 2 X and dropping the terms of second order in A, we get: 

v(r,t) = (-c 1 (t) 2 + |0(t) 2 )r (72) 
= E ext , 1 (r,t)-E a e d xtjl (r,t) (73) 

that substituted in eq . fl6~4]) yields E x y cl (r, t) = 0. 

"Strong" correlation limit The equilibrium density reduces to a 5-shell po — 
2/ {itro)5(r — ro/2) and we can treat the system classically considering the equation of motion 
of the separation Tyi between two classical point charges (see appendix C). Under the influ- 
ence of the external potential v exty i = ui{t) 2 r\ 2 /2 the equilibrium separation ri 2 oscillates 
according to the classical equation of motion 

fi2 = E exti i(t) - 3uq (n 2 - r ) (74) 

with E exti i(t) = —uf(t)ro and tq is the equilibrium separation in the absence of the external 
field (we use the linear response approximation). We can now define the "adiabatic external" 
field as the one that produces the same deviation from equilibrium as E ext ,i(t), under static 
conditions (f 12 = 0). That means that 

E&,i(t) = 3^ Q 2 (r 12 (t) - r ) (75) 
where r 12 (t) is the solution of (174]) . From this we can deduce that 

v, = E exttl {t) - Ef xtjl {t) (76) 

where v\ is the exact velocity field in this limit. This implies that E^ xt x (t) =0 in this limit. 
Non extreme cases 
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In this case the problem is to find a good approximation for E^ t (r, t). In the case we 
are considering E ext (r, t) = —u){t) 2 r so in order to have a simple and separable form for 
E"^ t (r, t), we can choose: 



very good results and p (r;i) is indistinguishable from p(r,t) within the numerical error. 
The results for the "high correlation" parameter uoq ~ 0.02 are less good. They can be 
improved using a(t) = uJq(1 + ecos(Qt)) and tuning the parameter e. In every case, also in 
these intermediate cases E^(r, t) ~ within the numerical error. 

We conclude that for this particular system the dynamical part of E xc (r, t) is almost 
negligible. However we caution that this is at least partly a special feature of the harmonic 
system studied here (see discussion in the following section) and should not be uncritically 
generalized to other systems. 



The comparisons performed in this paper between the exact xc potential of a two-electron 
harmonic atom and several approximate expressions for this quantity, constitute an ex- 
tremely severe test of the approximations in question. Aside from the exchange-only OEP, 
all the approximations considered are based on the homogeneous electron gas, and, therefore, 
are expected to be valid only for systems whose density is slowly varying on the scale of the 
local average inter-electron distance. This condition is certainly not satisfied by our model 
system - not in the weak correlation regime, in which the length scale of density variation 
coincides with the average inter-electron distance, and much less in the strong correlation 
regime, in which the latter greatly exceeds the former. In this light, the fact that the ALDA 
and GK produce potentials reasonably close to the exact ones, although qualitatively incor- 
rect at large distance from the center, should be regarded as an unexpected success of these 
approximations. 

Another surprising result of our study comes from the separation of E xc (r, t) into an 
"adiabatic" and a purely "dynamical" part. The somewhat counterintuitive result is that, 
in the case of a time-dependent harmonic external potential, the dynamical part of E xc 
is zero in the limits of weak and strong correlation and almost negligible in between. This 
happens at frequencies well above the first excitation threshold, where the density response is 
far from adiabatic. In the weak correlation regime, this result depends crucially on the form 
of the wavefunction in a parabolic potential. Therefore we don't expect the conclusion to be 
generalizable to other potentials. In the strong correlation regime, however, the reduction 
of the dynamics to harmonic oscillations about a classical equilibrium configuration appears 
to be a generic feature of the many-electron system. It is this feature that leads to the 
vanishing of E^(r, t) in this regime. 

These results throw some light on the surprising ability of the ALDA to give good results 
even outside its natural domain of validity (low frequency regime): in a system in which 




(</>) 2 r gives 



(77) 
(78) 



VI. DISCUSSION AND SUMMARY 
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the non-adiabatic corrections are small, a static functional of the density (such as the LDA 
xc potential) which works well in the static regime, is expected to give a reasonable time- 
dependent potential upon replacement of the static density with the time-dependent one. 

The recently introduced VUC approximation, contains a "dynamical" correction to 
ALDA (see eq. (|58|)) and, in the light of the exact behavior of E^(r, t) just underlined, 
it is interesting to notice that the "dynamical" part of VUC is, for this system, small, 
becoming almost negligible for strong correlation. 

In summary, we have found that, for this particular system, the "dynamical" part of E xc 
is almost negligible, and the ALDA, GK and VUC approximations work reasonably well 
at all coupling strengths (although the VUC underestimates E xc (r,t) for weak correlation). 
The OEP, as expected, is reasonable only for weak correlation. The main discrepancies are 
found to occur at small r and at large r (except for the OEP that has the exact asymptotic 
behavior). The question of whether these results are generalizable to more complex systems 
remains open. 

We acknowledge support from NSF Grant No. DMR-9706788 and from Research Board 
Grant RB 96-071 from the University of Missouri-Columbia. We thank C. A. Ullrich, S. 
Conti, R. Nifosi and C. Filippi for useful discussions. 

APPENDIX A: 

Making the change of variable R = 21Z, eq. (0) can be rewritten as: 

(-V 2 R + ^ 2 (t)R 2 )V CM (R,t) = iJU CM (R,t). (Al) 

Separating angular and radial coordinates as in eq. (|TTD , we obtain the radial equation: 
0^ 1 d \ tti^ 



v OR 2 RdR 4 jpw'K ' / dr 

Inserting into eq. ( |A2| ) the guess 

X n, m (R, *) = A(t)R m exp(B(t)R 2 )L r :(C(t)R 2 ) (A3) 

(a generalization of the corresponding static solution), we obtain the following equations for 
the time dependent coefficients: 

iC + 8BC + AC 2 = (A4) 
i- + AB + AmB - AnC = (A5) 
%B + AB 2 - ^o 2 {t) = 0. (A6) 



With the ansatz B = (i/A)(X/X), eq. flA~6|) becomes: 

X = -u 2 {t)X, (A7) 

the classical equation of motion for an harmonic oscillator. The solution X(t) can be written 

as 
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X(t) = \X(t)\e 1 ^ =Xx + tXz. 



(A8) 



Since Xn, m {R,t) must not diverge as R — > oo, we must impose that the real part of B(t) 
-Bsr + iB$. be negative. Using (|A8|) , B(t) can be written as: 



(A9) 
(A10) 



. , 1 dip i din \X\ 

Al) ~ ~A~dl + 4 It 

W idln\X\ 

~4|X| 2 + 4 dt 



where W = X^Xsr — X^X^ is a constant, being the Wronskian of two solutions of eq. 
In order to have a normalizable wave function not identically zero, we have then to impose 
that W > or equivalently that dcp/dt > 0. 



Requiring that C G 3? , from the real part of eq. ([A4]) and from eq. (|A10|) we get 

C(t) = (1/2)(W7|X| 2 ) (All) 



which also satisfies the imaginary part of eq . (|A4|) . Now we can solve eq. (|A5|) from which, 
integrating, we get: 



A(t) = A(0) ■ exp{-(m + 1) In 



X(0) 



i(2n + rn+l)((f>(t)-<f)(0))}. 



(A12) 



A(0) is determined by the normalization condition Jq°° \xnm 

(R,t)\ 2 RdR 



A(0) 



nl 



2 m (n + m)\ v dt 



,a<p \ 
\dt t=o) 



(m+l) 



(A13) 



Finally inserting all the expressions for the coefficients in eq. (|A3|), after some algebra we 
get: 



Xn,m 

(R,t;X) 



711 



V 2 m (n + m)\\dtJ 
i?-expf(-^ + 



/d(b\ m+1 / 

l-^) 2 exp(i(2n + m + l)(0(O)-0(t)) 



dln\X\ ^ R 2 ^ jm ( 
dt 



r>2 

4 ) n ^dt 2 



(A14) 



APPENDIX B: 

The solution of static radial equation for the RM channel can be written as: 



;bi) 



with 
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oo 

1 „2 



u(p) = e-*> ff^avpf, (B2) 

v=0 

p = ^fuT r r (B3) 
u^f (B4) 

s = - + Vl. (B5) 

The coefficients of the sum in ( |B2|) are related by: 

( a ( a + l)-/+I)o 1 = -^= (B6) 
4 ,AlV 



((u + s){v + s-l)-l + ~K - + - 1 - 2(1/ - 2 + s)]a„_ 2 = (B7) 

4 -./OV 



where e r is the part of the energy coming from the RM channel and ao is fixed by the 
normalization condition. Imposing the conditions a n _i ^ 0, a n = 0, a n +i = 0, the sum 
in (B2) can be made finite and the coefficients a v calculated. From these conditions we 
also obtain an expression for the energy e r = u r (2n + 2s — 1) and an expression (less 
straightforward) for u r . 

In our calculations we have truncated the sum in eq. ( |B2|) at n = 2 obtaining Uq — 1 
(weak correlation case) and at n = 5 obtaining cu = (25 — 3v / 33)/328 ~ 0.02 (strong 
correlation case). 



APPENDIX C: 

In the limit of strong correlation (eq. (|30D) and linear response regime the potential 
energy can be expanded up to second order about the classical solution and we can also 
approximate the momentum p with the radial component p r = —i(d/dr)f since, in this 
limit, the dynamic of the problem is basically confined in the f direction. The Hamiltonian 
of the relative motion problem eq. (||) can then be approximated as: 

H eff (t) = ^ + ^ 2 (r - r ) 2 - pE extjl (t)(r - r ) (CI) 

where p = 1/2 is the reduced mass, r = (1/ ' puj^) 1 ^ is the classical separation between 
electrons in the linear regime, oj 2 {t) = 3uq + oj\{t) and E ext ^{t) = —u>f(t)r can be viewed 
as an "external force". If we define r% = r — r Q , the deviation from the classical equilibrium 
position, we can use the change of variable r\ — y — Xo(t) so that the Hamiltonian becomes 

9 

H(t) = + 2& 2 y 2 - pu 2 x (t)y - px (t)y - pE ext>1 (t)y (C2) 

where we have dropped the irrelevant terms depending on the time alone. If we impose that 

x = -Cj 2 x - E ext>1 (t) (C3) 
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we get 



m) = f^^v (04) 

and the problem reduces to a 1-dimensional harmonic oscillator and can be solved exactly 
in a way similar to the one shown for the 2-dimensional harmonic oscillator in Appendix A. 
The general solution takes the form: 

*(y,t) = (-^)i(^)^i^ 2 e ^^)W°)-^) J ff n ((^)l 2/ ) (C5) 
2 n nU/TT at at 



where H n (x) are the Hermite polynomials and X(t) is a complex solution of the classical 
equation of motion 

X = -u 2 (t)X, (C6) 

= |X(t)|e^ (t) (C7) 

with a phase satisfying the condition d<j)/dt > 0. The solution of the original problem 
eq. ( |C1[ ) with n = is therefore 

#(r,t) = i r (^)3 e 5^ (0) -* ( * )) e"^ io(i){r " ro) e^$ (r " ro+Io( * ))2 e^o d *'(^ 2l o^ (C8) 

We stress that in the regime of high correlation A/ro — * 0, where A oc cJq the width of 
the Gaussian entering the solution (|U8|) , the wave function is concentrated around r (that 
justifies the approximation ( |C1| )) and tends to a 5-function in the extreme limit. 
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TABLES 



TABLE I. Parameters used in our numerical calculations. 







n 




A 


high corr. 


(25 - 3^/33)/328 « 0.02 


0.1 


» 4.2 


0.1 


low corr. 


1 


3.2 


3.2 


0.1 
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FIGURES 



FIG. 1. Lowest Floquet state electronic density for the weak correlation case (u>o = 1). The 
solid line is the exact static result while each of the broken lines corresponds to different times. In 
the inset we show the corresponding velocity field v(r,t). Each solid line corresponds to different 
times. 

FIG. 2. Same as Fig. 1 but for the strong correlation case (ojq ps 0.02). 

FIG. 3. Exact xc field E xc (r,t) for the weak correlation case (ojq = 1). The solid line is the 
static limit while each of the broken lines corresponds to E xc (r, t) at different times. Asymptotically 
Exc{r,t) f» — 1/r 2 , as in the static case. For comparison the static LDA result is also plotted. 

FIG. 4. Same as Fig. 3 but for the strong correlation case (ujq ~ 0.02). 

FIG. 5. Comparison between the exact first Fourier component E xc (r, Q) of the xc field and 
some of its most used approximations (weak correlation case, loo = 1). 

FIG. 6. Same as Fig. 5 but for the strong correlation case {ujq ps 0.02). 
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